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Abstract 

Summation-by-parts (SBP) operators are finite-difference operators that mimic integration by 
parts. This property can be useful in constructing energy-stable discretizations of partial dif- 
ferential equations. SBP operators are defined by a weight matrix and a difference operator, with 
the latter designed to approximate d/dx to a specified order of accuracy. The accuracy of the 
weight matrix as a quadrature rule is not explicitly part of the SBP definition. We show that SBP 
weight matrices are related to trapezoid rules with end corrections whose accuracy matches the 
corresponding difference operator at internal nodes. The accuracy of SBP quadrature extends to 
curvilinear domains provided the Jacobian is approximated with the same SBP operator used for 
the quadrature. This quadrature has significant implications for SBP-based discretizations; for ex- 
ample, the discrete norm accurately approximates the I? norm for functions, and multi-dimensional 
SBP discretizations accurately mimic the divergence theorem. 

Keywords: Summation-By-Parts operators, high-order quadrature, Euler-Maclaurin formula, 
end corrections, Gregory rules 



1. Introduction 

Partial differential equations (PDEs) are often solved numerically in order to approximate a 
functional that depends on the solution; for example, when computational fluid dynamics is used 
to estimate the lift and drag on an aerodynamic body. For integral functionals, such as lift and 
drag, a quadrature rule may be needed to numerically integrate the discrete solution. When we 
are free to choose the quadrature weights and abscissas, Guassian quadrature is often the optimal 
choice. However, the choice of quadrature rule is less clear for the uniform grids that arise in 
finite-difference methods. 

In this paper, we investigate a quadrature rule that is particularly well suited for high-order 
summation-by-parts (SBP) finite-difference methods SBP operators lead to linearly time-stable 
discretizations of well-posed PDEs, and they have been used to construct efficient discretizations 
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of the Euler 0, 0], Navier-Stokes 0-0], and Einstein equations 0|. The high-order quadrature 
in question is based on the weight matrix that forms part of the definition of SBP operators. 
This result is somewhat surprising, because the accuracy of the quadrature induced by the weight 
matrix is not explicitly part of the SBP definition. To our knowledge, the relationship between 
SBP operators and quadrature has not been discussed previously in the literature. 

In the context of high-order finite-difference methods, including those based on SBP operators, 
several classical quadrature rules are available to accurately evaluate integral functionals; for exam- 
ple, composite Newton-Cotes rules and Gregory-type formulae. Why introduce a new quadrature 
rule based on SBP weight matrices? While accuracy is important, we may also want the func- 
tional estimate to obey some property or properties of the true functional, and this is the value of 
SBP-based quadrature. 

Consider the volume integral of the divergence of a vector field over a compact domain. The 
resulting functional is equivalent to the flux of the vector field over the domain's boundary, in light 
of the divergence theorem. This is a fundamental property of the functional that we may want a 
discretization and quadrature to preserve. We say a functional estimate respects, or mimics, the 
divergence theorem if 1) it is accurate, and 2) the discrete quadrature over the volume produces a 
discrete quadrature over the surface. 

In general, classical quadrature rules for uniformly spaced data will not mimic the divergence 
theorem in the above sense when applied to an arbitrary high-order finite-difference approximation 
of the divergence; typically, they will satisfy the first but not the second property. In contrast, 
we will show that an SBP discretization does mimic the divergence theorem when numerically 
integrated using its corresponding weight matrix. 

The paper is organized as follows. Section [5] introduces notation and formally defines SBP 
operators. Section [3] presents the main theoretical results. In particular, we derive conditions on 
the quadrature weights for the class of trapezoid rules with end corrections. These conditions are 
used to establish the accuracy of SBP-based quadrature. Subsequently, we consider the impact 
of coordinate transformations on SBP quadrature and show that the quadrature remains accurate 
on curvilinear multi-dimensional domains. In Section [3] we verify the theoretical results with 
several numerical examples. The implications of SBP quadrature are summarized and discussed in 
Section [5j 

2. Notation and definitions 

We try to remain consistent with the notation used by Kreiss and Scherer in their original work 
[l[, as well as Strand's subsequent work [§]. 

The interval [0, 1] is partitioned into n + 1 evenly spaced points x v = vh, v = 0, 1, . . . , n, with 
mesh spacing h = 1/n. Finite intervals other than [0,1], as well as nonuniform node spacing, 
can be accommodated by introducing an appropriate mapping (see Section \'6.2\) . For arbitrary 
U(x) G C p [0, 1], we use u v = U{x v ) to denote the restriction of U to the grid x v . 

Definition 1 (Summation-By- Parts Operator). The matrix D G R( n+1 ) x ( n+1 ) is a summation- 
by-parts operator for the first derivative on the mesh {a^}™^ if it has the form 

D = H~ l Q, 

where the weight matrix H G ]j( n + 1 ) x ( n + 1 ) ig a symmetric-positive-definite matrix, and Q G 
M(n+ i )x{n+ i) satigfies 

Q + Q T = diag(-l,0,0,...,l). 
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Furthermore, the truncation error of the difference operator D in approximating d/dx is order h 2s 
at the internal nodes, {x v }^Z^, and order h T at the boundary nodes, {x„}^~q and {x v }™ =n _ r+1 , 
where r, r, s > 1 . 

In other words, the SBP operator D approximates d/dx and has a particular structure. In 
general, the order of accuracy of the difference stencil at internal nodes is different than the order 
of accuracy of the stencil at boundary nodes. The even order of accuracy 2s for the internal nodes 
is a consequence of using centered-difference schemes, which provide the lowest error for a given 
stencil size. For a 2s-order accurate scheme, the derivative at the internal nodes is approximated 
as 

xl^) ~ ^1 - Uw-v), r <w <n-r, 



where the coefficients are defined by (see [9j, for example) 

(-ir+w 



a:. 



v(s + v)l(s — v)\ 



The following lemma from [8| lists some identities that the a v satisfy; these identities will be useful 
in our subsequent analysis. 

Lemma 1. The coefficients a v that define a 2s-order accurate SBP operator at internal nodes 
satisfy 



v=l 



0, j = l,2,...,s-l. 



We turn our attention to the weight matrix H, which is the focus of this paper. Since H 
is symmetric-positive-definite, we can use it to define an inner product and corresponding norm 
for vectors. Let u, z € M. n+1 be two discrete functions on the grid nodes, i.e. u v = U(x v ) and 
z v = 2(x v ). Then 

(u,z)h = u T Hz, and \\u\\% = (u, u)h, 

define the H inner product and H norm, respectively. Using the SBP-operator definition and the 
H inner product, we have 

(u,Dz)h = -(Du,z) H -u z + u n z n . (1) 
Equation (pQ) expresses the fundamental property of SBP operators and is the discrete analog of 



This property of SBP operators is what leads to energy-stable discretizations of partial differential 
equations. However, while (pQ) is analogous to integration by parts, it remains to be shown that 
([!]) is an accurate discretization of ([2|). 
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In this work, we will consider H matrices with the block structure 

(H L \ 

H = h I , (3) 




where Hl,Hr £ M rxr are symmetric-positive-definite matrices. Assuming that ii/, and Hr are 
dense matrices — the so-called full-norm case — Kreiss and Scherer [l[ established the existence of 
SBP operators that achieve an order of accuracy of r = 2s — 1 at the boundary with r = 2s. Strand 
0] showed that 2s — 1 accuracy can be maintained at the boundary in the case of a restricted-full 
norm, which uses 

with H L , H R € Rfr-iM'- 1 ) and r = 2s + 1. 

In general, SBP weight matrices of the form ([3]) satisfy the compatibility conditions described 
in the following proposition [l|; these conditions will be used later to establish the accuracy of 
quadrature rules based on full and restricted-full H matrices. 

Proposition 1. Let H € R( n+1 ) x ( n+1 ) be an SBP weight matrix with the block structure ([3]). Then 
Hl satisfies 

jejH Lej -i + ie]H L ei-i = -{-r) i+j + J iJ; < i, j < r, 
where ej = (— 1) J (r 5 ' (r — 1) J • • • with the convention e_i = 0, and 

-v-i 



Ji,j — Ct v 



v=l 



it^(ti; — v) 4 + w % {w — v) 3 



.to=0 



i + i > i. 



Kreiss and Scherer also showed that it is possible to define SBP operators with diagonal H 
matrices, i.e. 

H L = diag (A , Ai, . . . , A r _i) 
H R = diag (A r _i, . . . , Ai, A ) 

with Aj > 0. These "diagonal norms" are important because, unlike full and restricted-full norms, 
they lead to provably stable PDE discretizations on curvilinear grids However, diagonal-norm 
SBP operators are limited to r = s accuracy at the boundary when the internal accuracy is 2s. 
Consequently, the solution accuracy of hyperbolic systems discretized with such SBP operators is 
limited to order s + 1 [lit]. Nevertheless, one can show that functionals based on the solution of 
dual consistent diagonal-norm SBP discretizations are 2s-order accurate H, 13]. 

When the weight matrix H is diagonal, Kreiss and Scherer Q showed that its elements are 
defined by the relations in following proposition. 

Proposition 2. Let H € R( n+1 ) x ( n+1 ) be a diagonal SBP weight matrix with r = 2r = 2s. Then 
the diagonal elements X v of Hi and Hr satisfy the relations 

v -i / (rP-(-l)% j = l,2,...,2s-l 

3 l^ X v{r -vy - 



v=0 



where (3j is the j th Bernoulli number. 
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3. Theory 



3.1. One- dimensional SBP Quadrature 

To establish the accuracy of SBP-based quadratures, we need the following theorem that places 
constraints on the coefficients of a certain class of quadrature rules for uniformly spaced data; 
specifically, the trapezoid rule with end corrections. The theorem is a direct consequence of sub- 
stituting finite-difference approximations into the Euler-Maclaurin sum formula. 

Theorem 1. Consider a set of n + 1 uniformly spaced points, x v = vh, v = 0,1,..., n, with 
constant mesh spacing h = 1/n. A quadrature of the form 



'r— 1 n—r i — 1 



Z(u) = h I } j a v u v + u v + cr v 



w=0 v=r v=0 



IS 



a q- order accurate approximation of f^Udx forlA € C 2m+2 [0, 1], where q — 1 < r and q < 2m+2, 
if and only if the coefficients {o~ v } r v ZJo satisfy 

r-1 

jY,a v (r-vy- 1 = ri-(-lY0 j) j = 1, 2, . . . , q - 1, (4) 

v=0 

Proof. Consider the Euler-Maclaurin sum formula applied to U(x) [3]: 

/ U{x)dx = hY^u v + Y^^h k (ut 1) -{-l) k u^t 1) )+E 2mi (5) 
■'° u=0 fc=l ' 

where it!* 1 ^ = D^ k ^ 1 ^U(x v ), 2m < q < 2m + 2, and the error term is given by 

a u2m+2 

^ 2m - (2m + 2)! D 

with £ G (0,1). Suppose the function derivatives at x = and x = 1 are replaced with finite- 
difference approximations involving the first r and last r internal points, respectively. Moreover, 
assume that the approximation to ui k is accurate to 0(h q ~ k ), where q — 1 < r; consequently, 
the approximations are exact for polynomials up to at least degree q—1. Let {di k denote 
the coefficients defining the finite-difference approximation of , such that 

r-1 r(fc-l) 

Substituting the finite-difference approximations into ©, and noting that the coefficients for odd 
derivatives must be negated at x = 1, we find 

1 n 2m a r— 1 ^(fc— 1) 

" ' 2m+2\ 



/ u(x) dx = hY,u v + J2j[ hk Yl ijrr ^ + n «-) + + 

i>=0 fc=l ' ?;=0 

/ 1 — 1 n—r r— 1 \ 

= /t I ^ + ^ U v + ^ VvUn-v + 0(h q ) 



\v=0 v=r v=0 
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where 

2m 

^ = 1 + Et^ _1) ' v = 0,l,...,r-l. (6) 

k=i 

Next, we will show that these a v are the same ones that satisfy (HJ), a set of q — 1 conditions that 

(k) 

are independent of the 5{, . Substituting the above expression for a v into (jU), we find 

r— 1 r— 1 2m „ r— 1 

v=0 v=0 k=l ' v=0 

j = l,2,...,q-l. (7) 
The first term on the right-hand side can be recast using the sum of powers formulsH: 

3 E( r - v ) 3 ~ l = rj + E(- X ) fc u ^~ k - ( § ) 

v=0 k=l ^ ' 



For the second term, we recognize that (r — v is the discrete representation of the polynomial 
Pj-i(x) = h~^~^'{rh — x)i~ l ] therefore, since the finite-difference approximations are exact for 
polynomials of degree q — 1, we have 



r-1 



u=0 

and 



0, k > j, 



Zm „ r-1 3 / \ 

iE§E e-^ (»■ - 1 = E(- 1 ) fc - 1 U j o) 

fc=l ' v=0 k=l ^ ' 

Substituting (jSJ) and ([9]) into (J7]), and recalling that the odd Bernoulli numbers greater than one 
are zero, we have 

j ]T a v (r - vy~ l = r j + £(-1)* (£) (3 k r^ k + ^(-l)^ 1 (£) ftr*"* 

v=0 k=l ^ ' fc=l ^ ' 

for j = 1, 2, . . . , q — 1. Thus, we have shown that the a v satisfy when the quadrature is g-order 
accurate. 

We need the general solution of @ to show that these conditions are sufficient for the quadra- 
ture to be g-order accurate. We have already shown that ^ is a particular solution of the linear 
equations (UJ), so we need to determine the form of the homogeneous solution, i.e. the null space 
of the matrix on the left side of (dJ . 

As noted above, (r — vy~ l is simply the polynomial pj—\(x) = h^^^ l \rh — x)i~ l evaluated at 
the nodes. The derivative operator Z)^' -1 ) with q < k < r will annihilate pj-i(x), since j < q — 1; 
therefore, any finite difference approximation that is a consistent approximation of h k ~ 1 D( k ~ 1 \ 



3 We use the sum of powers formula that is consistent with /3i = — | 
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q < k < r, will annihilate pj-\{x v ) = (r — v)i 1 . If we let {fii k }JJ = q 1 denote the coefficients of 
such a finite difference approximation, then the general solution to (jlj) can be written as 

^i+E^ + E^-/- 11 - (io) 

fc=l fc=g 

where {70,71, • • • , 7r-<y} parameterizes the null space. When r = q — 1, the null space is trivial, 
and the second sum does not appear in (fTUj) , 

Substituting the general solution into the quadrature yields 



'r— 1 n— r r— 1 



l(u) = h I cj^u„ + u„ + 



\i>=0 i>=r n=0 / 

n 2m „ r— 1 ^(fc— 1) 

to X! Uv + X 77^ X -jj^t K + 

u=0 k=l ' v=0 

r r— 1 

fc=q v=0 

f 1 U{x) dx + 0(h q ) + V 7 fc- g to fc (u?^ + (-l) fc 4 fc_1) 
Jo k= q V 

1 

W(x)dx + 0(/i 9 ). 



Therefore, we have shown that is sufficient for the quadrature to be g-order accurate, which 
completes the proof. 

If we choose q — 1 = r, Theorem [1] provides a closed set of equations for constructing high- 
order quadrature rules for uniformly spaced data with equal weights on the internal points. More 
generally, we may choose q — 1 < r, in which case the additional degrees of freedom can be used to 
achieve other objectives. For example, setting ao to zero, so that only strictly internal points are 
used. 

Theorem Q] encompasses many existing quadrature rules, including the Gregory class of for- 
mulae, and it could be used to construct an unlimited number of novel trapezoid rules with end 
corrections. However, our interest in Theorem [1] is not in constructing new quadrature rules, but 
in its consequences for SBP weight matrices. 

Corollary 1. Let H be a full, restricted- full, or diagonal weight matrix from an SBP first-derivative 
operator D = (H~ 1 Q), which is a 2s -order- accurate approximation to d/dx in the interior. Then 
the H matrix constitutes a 2s -order- accurate quadrature for integrands U £ C 2s . 

Proof. For diagonal SBP weight matrices the result follows immediately from Proposition^ since 
flU), with q = 2s, is a subset of the equations that define the X v . For the full and restricted- full 
weight matrices, consider the relations in Proposition [1] with j < r = 2s — 1 and i = 0: 

r— 1 r— 1 s v—1 

j'EE K w {-iy-\r - wy- 1 = -(- r y + £ y + ( w - v y] . 

v=0 w=0 v=l w=0 
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Multiplying the left and right sides by (— 1) J 1 , using the symmetry of the h vw , and swapping 
summation indices on the left side, we find 



r-l 



v-l 



j ^2<r v (r- vy- 1 = r j + 1 ^ a v ^ [w j + 



w — V) 



v=0 



v=l 



io=0 



where a v is identified with Ylw=o^vw The secon d term on the right-hand side can be simplified 
using the accuracy conditions of the a v (Lemma [T]) and the formula for the sum of powers. 



v-l 



v=l 



w=0 



v=l 

s 



W J + 



w=l 



v=l 



v + " + ± 

to=0 



10 = 1 

J + l 



3 + 1 



w 



j+l-w 



Thus, we have 



-2. 3 = 1 

0, j = 3, 5, . . . , r, 

[ pj, j = 2,4,...,r-l 
Pi- 



r-l 

j ~ v)^ 1 = r j - (-lyfr, l<j<r, 

v=0 



and Theorem [T] implies that full and restricted-full SBP weight matrices are quadrature rules 
accurate to r + 1 = 2s. 

3.2. SBP Quadrature and Coordinate Transformations 

Curvilinear coordinate systems are often necessary when solving PDEs on complex domains. 
Like most finite-difference schemes, SBP operators are not applied directly to the nodes in physical 
space. Instead, a coordinate transformation is used to map points in the physical domain to 
points on a Cartesian grid, and the SBP operators are applied in this uniform computational 
space. However, this coordinate transformation introduces geometric terms whose impact on the 
accuracy of the quadrature rule is not clear. 

We begin by considering the one-dimensional case. Let T(x) = £(x) be an invertible transfor- 
mation of class C 2s that maps Q x = [a,b] to = [0, 1]. For Li G L 2 (Q X ), the change of variable 
theorem implies 



U dx 



UJdi, 



;n) 



where J = ^| is the Jacobian of T _1 . 

We are interested in the accuracy of SBP quadrature in the computational domain, so we 
consider the discrete equivalent of the right-hand side of (jlip . In general the mapping will not 
be explicitly available, so the Jacobian must be approximated. As we shall see, to retain the 2s- 
order accuracy of SBP quadrature, it is critical that the derivative that appears in the Jacobian 
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be approximated by the same SBP difference operator that defines the norm. Thus, if x € M Tl+1 
denotes the coordinates of the nodes in physical space, the SBP approximation of (jlip is given by 



u 1 HDx = u 1 Qx. (12) 

The following theorem confirms that this discrete product is a 2s-order accurate approximation of 
the integral (fTT]). 



Theorem 2. Let D = H~ 1 Q be an SBP first derivative operator. Then 

(z, Du)h = z T Qu 
is a 2s -order- accurate approximation to the integral 

[ Z—dx, 
Jo ax 

where Z<§ G C 2s [0,l}. 

Proof. Using SBP-norm quadrature we have 

£ Z^-dx = (z,u') H + 0(h 2s ), 

where u' denotes the analytical derivative dU/dx evaluated at the grid nodes. The result will follow 
if we can show that 

(z,u') H = (z,Du) H + 0(h 2s ). (13) 

The expression on the left of (|13[) is simply a quadrature for the integrand y = Z-jj^lA. Con- 
sequently, it is sufficient to show (fT3j) is exact for polynomial integrands of degree less than 2s. 
Let 

w i = [^0 x l ' ' ' X n] 

be the restriction of the monomial x l to the grid. We will consider 

z = Wi, u = Wj, and u = jvjj-i, 

with i+j < 2s. 

First, suppose j < s. In this case, an SBP operator (including those with diagonal-norms) is 
exact for Wj giving 

Du = Dujj = jwj-i = u', 

and substitution into (fT3|) yields (z,u')h = (z,Du)h- 

Next, to show that (|13|) is exact for j > s, the roles of z and u will be reversed. Here, since 
j + i < 2s, we must have i < s, and the SBP operator becomes exact for Wi\ 

Dz = Dwi = iwi-i = z . 
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Using this exact derivative and the properties of SBP operators we find 



{z,Du) H = z T H{H- 1 Q)u 

= z n u n - z u - z T Q T u 




(u,Dz) H 
{u,z') h 



= f 4- (UZ) dx - [ 
Jo dx Jo 





Thus we have shown that the expression (z, Du)h is also equal to the exact integral when j > s 
and i + j < 2s. This completes the proof. 

For multidimensional problems on curvilinear tensor-product domains, SBP operators are ob- 
tained from the one-dimensional operators using Kronecker products. To extend SBP quadrature 
to these domains, we need only apply Theorem [2] iteratively over the individual coordinate direc- 
tions. We provide a sketch of the proof here and direct the interested reader to [131 ] for the details 
of the two-dimensional case. Consider the change of variable theorem in d dimensions: 



where J is the Jacobian of the mapping (more precisely, the determinant of the Jacobian). As in 
the one-dimensional case, the mapping and integrand must be sufficiently differentiable (class C 2s ) 
for the quadrature to remain 2s-order accurate. An important observation is that the Jacobian 
consists of a sum of terms of the form 



in which none of the indices i,j,...,k are equal. Because the indices of the computational coor- 
dinates are also distinct, Theorem [2] can be applied one dimension at a time (i.e., as an iterated 
integral). For example, we can consider dimension and apply Theorem [2] to the integral 



corresponds with Z. Repeating this process over the remaining coordinate directions and terms in 
the Jacobian yields the desired result. 




dxi dxj dxk 




where Xi corresponds with hi in the theorem, and 
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3.3. SBP Operators and the Divergence Theorem 

Using the above results, one can show that SBP operators mimic the d-dimensional divergence 
theorem to order h 2s on curvilinear domains that are diffeomorphic to the ci-cube. We will consider 
the two-dimensional case; the extension to higher dimensions is straightforward. 

In two-dimensions, the divergence theorem is 



dT j_ 96 A A 

— + — dxdy 
ox oy 



[Tdy + Qdx) 



(14) 



where d$l x is the piecewise-smooth boundary of £l x , oriented counter-clockwise. Applying the 
coordinate transformation, we find 



dT j. 95 A A 

+ — dxdy 
oy 



dx 



dT 

dx 



dg 

dy 



J dxdy 



II 



dT u_ " A, A 



(15) 



where we have used the metric relations 15j, [lfj to obtain the components 



J 
G = J 



dx dy 
\ox dy 



dy_ 
drj 
dy 
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dx 
— ( 
dr]' 



, dx 

di T+ M G - 



(16) 
(17) 



In light of (|15|) . we need only show that SBP discretizations obey the divergence theorem to order 
h 2s in the simpler computational space: 



dT JU dG A, A 

-W + d^ d ^ dr] 







(18) 



The reader may object to this simplification, since T and Q contain derivatives that depend on 
the geometry and must be approximated. However, if the partial derivatives of x and y appearing in 
(fT6|) and (fTT|) are approximated using the same SBP operators as found in the discrete divergence 
theorem, then Theorem [2] can be applied. This follows because the same difference operator is 
never applied twice in the same coordinate direction (e.g., d/d£ is applied to J 7 , which contains 
only partial derivatives with respect 77). 

For simplicity, assume that the square fi^ is discretized using n + 1 nodes in both the £ and rj 
directions. Thus, the nodal coordinates are given by 



0<j,k< n}. 



If the nodes are ordered first by j and then by k, one-dimensional SBP operators can be used to 
construct the two-dimensional difference operators 



D C = (I ® D) 



and D v = (D®I), 
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where (g> denotes the Kronecker product, D = H~ 1 Q is the one-dimensional SBP operator, and I 
is the (n + 1) x (n + 1) identity matrix. Similarly, (H (g) H) defines the SBP quadrature for the 
two-dimensional set of points. Let B = diag (—1,0, 0, . . . , 1), so that we may write Q + Q T = B. 
Finally, let / and g denote the restriction of the functions T and Q, respectively, to the grid points, 
and let c = [l 1 • • • l] denote the constant function 1 restricted to the grid. 

With the two-dimensional SBP operators suitably defined, we can discretize the left-hand side 
of (USD: 



c T {H®H) (I®D)f + {D®I)g 



c T (H®Q)f + c T (Q®H)g 

c T (H®(B- Q T )) f + c T ((B - Q T ) ®H)g 

n n 

hii(fn,j — fo,j) + hii{9i,n — 9i,o), 
j=0 i=0 



(19) 



where we have used c T (Q T ® H) = c T (H <£> Q T ) = (constants are in the null space of D = H T Q). 
We highlight two significant facts regarding (|19|) . 

1. It is a 2s-order accurate approximation of the right-hand side of (|18l) . 

2. It depends only on the terms of / and g that fall on the boundary. 

Constructing a scheme that satisfies either one of these properties may not be difficult; however, 
few high-order schemes satisfy both 1 and 2 simultaneously. This is what we mean when we say 
the SBP operator mimics the divergence theorem. 



4. Examples 

4-1. One- dimensional Quadrature 

To illustrate the basic theory, we use the weight matrices from several common SBP operators 
to integrate a simple function. We consider three SBP operators with diagonal weight matrices 



and one SBP operator with a full norm. The diagonal operators are taken from Diener et al. [17] 
and are denoted by diag-r-2s, where r and 2s indicate the truncation error at the boundary and 
interior, respectively. The full norm operator can be found in [8| and is denoted full-r-2s. The 
boundary weights a v for all four operators are listed in Table [TJ for the diagonal norms o~ v — A^, 
whereas for the full norm a v = Y^w=o h vw . 
Consider the definite integral 

1= f U{x)dx 
Jo 

= / (4vr) 2 a;cos(47rx)dx (20) 
J o 

= — 4-7T COS (47r). 

To assess the accuracy of the SBP quadrature rules in Table [H we perform a grid refinement 
study based on the integral (f20j) and using n £ {16, 32, 64, 128, 256, 512}. Table [2] lists the rates of 
convergence for the quadrature rules. For n > 16, the rate of convergence is calculated from 

(21) 




Table 1: Boundary quadrature weights corresponding to some SBP weight matrices. 



SBP operator 


T 


2s 




or 


C"2 


"3 


<T 4 


C"5 


diag-1-2^ 


1 


2 


l 

2 












diag-2-4 


2 


4 


17 
18 


59 
48 


43 
48 


19 
48 






full-3-4 


3 


4 


43 
144 


67 
48 


35 
48 


155 
144 






diag-3-6 


3 


6 


13649 
43200 


12013 
8640 


2711 
1320 


5359 
4320 


7877 
8640 


43801 
43200 



the trapezoidal rule 



Table 2: Rates of convergence for the SBP quadrature rules in Table Q] applied to (ffil)) . 









n 






SBP operator 


32 


64 


128 


256 


512 


diag-1-2 


2.0113 


2.0028 


2.0007 


2.0002 


2.0000 


diag-2-4 


4.4978 


4.4148 


4.2182 


4.1019 


4.0473 


full-3-4 


4.1973 


2.9369 


3.7072 


3.8876 


3.9510 


diag-3-6 


5.7050 


6.8942 


6.9378 


6.7651 


6.5472 



where E n = X — c T Hu, with c T = (l 1 ... l) , is the error using n + 1 nodes. In all cases, the 
errors converge to zero at the expected asymptotic rate of 2s. 

Figure [T] plots the errors E n versus a normalized mesh spacing. This figure reminds us that 
schemes with the same order of accuracy can produce different absolute errors: the diag-2-4 opera- 
tor is almost an order of magnitude more accurate than the full-3-4 operator for n > 64. However, 
further analysis is required before we can characterize the relative performance of these schemes 
more generally. 

4-2. Multi- dimensional Quadrature on a Curvilinear Domain 

As shown in Section 13.21 SBP quadrature retains its theoretical accuracy on curvilinear do- 
mains provided the Jacobian of the transformation is approximated using the corresponding SBP 
difference operator. To verify this, we consider the domain 

Q x = {( x , y) G R 2 | 1 < xy < 3, 1 < x 2 - y 2 < 4}, 

and the integral 

-7- ll\ 2 2\ 1 -' x ' 2 +v' 2 ( X V — 1\ 

I = (x + y )e 3 sin ( — - — I dxay 

= 3(1 - e^ 1 )(l - cos(l)). (22) 
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♦■" full-3-4 
-A- - diag-3-6 



10 



»-1 



10" 



(n/16)" 

Figure 1: Errors of the SBP-based quadrature rules applied to (|20|l . 



To compute this integral numerically, we introduce a computational domain based on the coordi- 
nates 



x 



y 



l 



and 7] 



xy — 1 



For a given n £ {16, 32, 64, 128, 256, 512}, we divide £ and r/ uniformly into n + 1 points to produce 
a Cartesian grid on the square VL^ = [0, l] 2 . The physical coordinates x and y are evaluated at each 
computational coordinate, and these are used to compute the integrand in ([22]) . which we denote 
by /. The grid for n = 32 is shown in Figure [2j 

The Jacobian of the transformation is approximated using 



J=[{I® D)x] o [(D ® I)y] - [(I ® D)y] o [(D ® i>] , 



(23) 



where o denotes the Hadamard product (the entry- wise product, analogous to matrix addition). 
We have assumed that the nodes are ordered first by £ and then by r/ , so we can construct the 
two-dimensional derivative operators using Kronecker products of the one-dimensional operator D 
and identity matrix /. 

For a given n, the SBP-based approximation of ()22[) is given by 



l n = J T {H®H)f, 



(24) 



and the error in the quadrature is E n = X — X n . As before, the order of convergence for n > 16 
is estimated by q n given by (|2~T|) . Figure [3] plots E n and Table [3] lists q n for the diagonal-norm 
SBP operators listed in Table [TJ We focus on the diagonal-norm SBP operators here, because their 
accuracy is less obvious; the derivatives used to approximate the Jacobian are only s-order accurate 



14 



1.6 



1.4 



1.2 



y 1 



0.8 



0.6 



0.4 

1.2 1.4 1.6 1.8 2 2.2 2.4 
X 

Figure 2: Example grid for Q x with n = 32. 



Table 3: Rates of convergence for the diagonal- norm SBP operator approximation of (|22l) . 



n 



SBP operator 


32 


64 


128 


256 


512 


diag-1-2 


2.0911 


2.0453 


2.0226 


2.0113 


2.0056 


diag-2-4 


4.3283 


4.1583 


4.0768 


4.0374 


4.0093 


diag-3-6 


7.0799 


6.7941 


6.2253 


2.1274 


-0.7390 


mixed 


3.3170 


2.0521 


2.7215 


2.8863 


2.9484 



at the boundary. Nevertheless, as predicted by the theory, Table [3] shows that the quadrature for 
these schemes remains 2s-order accurate. Note that the errors for the diag-3-6 scheme are corrupted 
by round-off error for n = 256 and n = 512, which explains the suboptimal values of q n for these 
grids. 

We have also included results for a mixed scheme in Table [3] and Figure [3l This mixed scheme 
uses the diag-3-6 SBP operator to evaluate the derivatives in the Jacobian (f23|) and the diag-2-4 
operator to evaluate the quadrature ()24p . The results show that the mixed scheme has an asymp- 
totic convergence rate of only 3. Thus, despite a more accurate approximation of the Jacobian, the 
mixed scheme produces a less accurate X n than the scheme using the diag-2-4 operator for both 
the Jacobian and quadrature. This illustrates the importance of using the same operator to obtain 
the theoretical convergence rate. 

4-3. Discrete Divergence Theorem 

In this final example, we verify that SBP operators mimic the divergence theorem accurately. 
Specifically, we wish to show that when the divergence of a vector field is discretized using SBP 
operators and then integrated using the corresponding SBP quadrature rule, the result depends 
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Figure 3: Errors of the SBP-based quadrature rules in approximating Q22p . 



only on the nodes along the boundary and is a 2s-order approximation to the surface flux. 

We adopt the same domain £l x and coordinate transformation as in the previous example. A 
vector field (J 7 , Q) is defined by 

x ( l-xy \ f 2n(x 2 - y 2 2y ( xy-l \ 7 . (n{x 2 -y 2 -l) 

y)= 77 ex P —5 — cos o + -5- — ^— sm 



2 1 \ 2 j \ 3 y 3 V 2 / V 3 
«,.,) . -|exp (1^2) cos ( 2 ^7 2 - 1 » ) + f (^)' S in "»* " 

The analytical value of the divergence of (J 7 , £/) integrated over the domain fi x is 

The discrete divergence is evaluated in computational space using approximations for T and Q. 
In particular, the derivatives of the spatial coordinates that appear in (|16p and (|17p are approxi- 
mated using SBP operators. Therefore, at the nodes, J- and take on the values 

/ = [(D ® I)y] of-[{D® I)x] o g, 
g = - [(I ® D)y]o f + ® D)x]o g, 



where / and g denote the values of T and Q evaluated at the nodes. 
The SBP approximation of X is given by (see (|15p ) 



T 



(I ® D)f + (D ® I)g 
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Table 4: Rates of convergence for the diagonal- norm SBP operator approximation of an integrated divergence 
field. Round-off errors are contaminating the estimates for diag-3-6 with n = 256 and n = 512 

n 

SBP operator 32 64 128 256 512 

diag-1-2 2.0909 2.0453 2.0226 2.0113 2.0056 

diag-2-4 3.7201 3.7862 3.9000 3.9532 3.9758 
diag-3-6 7.5935 7.2371 7.8361 5.0507 -2.1760 



Table H] lists the estimated order of accuracy q n based on E n = I — X n for the three diagonal- 
norm SBP operators diag-1-2, diag-2-4, and diag-3-6. As predicted, the SBP discrete divergence 
integrated using (H ®H) is a 2s-order accurate approximation to I. Moreover, in light of ()19p . we 
know that I n depends only on the boundary nodes (This has been confirmed by calculating the 
right-hand side of (|19[) and showing that it equals I n to machine error). 

5. Discussion and Conclusions 

We have shown that the weight matrices of SBP finite-difference operators are related to trape- 
zoid rules with end corrections. We make no claim regarding the optimality of these SBP quadrature 
rules with respect to existing schemes on uniformly spaced grid points. However, the result has 
significant implications for SBP discretizations of PDEs, which we list below. 

• The SBP energy norm, which is frequently used in the stability analysis of SBP-based PDE 
discretizations, is a 0(h 2s ) accurate approximation of the I? norm for functions on [0, 1]. 

• The summation- by-parts property, equation ([1]), is a formal and accurate representation of 
integration by parts, equation @. More generally, multi-dimensional SBP discretizations 
using Kronecker products mimic the divergence theorem, i.e. the weight-matrix quadrature 
applied to the discrete divergence produces an accurate quadrature of the flux over the domain 
boundary in which no interior points are involved. 

• Diagonal-norm SBP operators have s order-accurate boundary closures when the interior 
scheme is 2s-order accurate. This limits numerical PDE solutions to s + 1 order accuracy 



111 ]: however, in 13|] we show that a dual consistent SBP discretization leads to super- 
convergent 2s-order-accurate functionals, if the corresponding SBP quadrature rule is used 
to calculate the functional (see also [T3] ) . 

In light of these observations, the SBP weight matrix appears to be the natural quadrature rule 
for evaluating functionals from corresponding SBP discretizations. 

Other than its asymptotic form, the leading error term for SBP-based quadrature remains 
unknown, so we cannot make general statements regarding the relative accuracy of different weight 
matrices. Although not pursued here, characterizing the leading error term should be possible by 
finding the finite-difference schemes corresponding to the Sy in Theorem [TJ 
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